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~^The  evolution  of  numerical  techniques  for  solving  problems  in 
Fluid  Dynamics  is  followed,  in  outline,  from  the  days  when  Digital 
Computers  were  first  available,  at  the  end  of  the  Second  World  War, 
to  the  present  time,  when  the  Computer  Aerodynamic  Simulator  is 
being  assembled.  In  this  period  the  range  of  numerical  methods 
has  been  broadened  five  fold,  while  the  speed  and  capacity  of 
computers  have  increased  by  several  orders  of  magnitude.  Two 
areas  close  to  the’  author's  interests  are  selected  to  illustrate 
these  changes.  The  first  concerns  the  extension  of  the  Method 
of  Integral  Relations  to  apply  to  laminar  and  turbulent  boundary 
layer  problems,  including  internal  flows,  separated  flows  and 
turbulent  mixing  flows.  The  second  area  deals  with  unsteady 
inviscid  compressible  flow  in  one  or  more  dimensions  and  a  dis¬ 
cussion  is  given  of  the  relative  merits  of  Godunov  and  Glimm 
techniques .  - 

r-N 

INTRODUCTION 


The  need  for  numerical,  as  distinct  from  analytical,  methods  to  solve  problems  in 
Fluid  Dynamics  first  arose  when  problems  in  Gas  Dynamics  presented  themselves. 
These  are  connected  with  combustion  and  explosive  processes  and  are  essentially 
non-linear.  The  eauations  of  motion  governing  such  problems  are  hyperbolic  ana, 
in  principle,  could  be  solved  numerically  by  the  Method  of  Characteristics.  Th’s 
has  a  long  history  going  back  to  the  pioneer  paper  by  Massau  (1900).  The  founda¬ 
tions  of  the  method  are  laid  out  in  Courant,  Friedrichs  and  Lewy  (1923)  and 
developed  for  application  in  Courant  and  Friedrichs  (1978). 

Two  difficulties  arise  in  applying  the  Method  of  Characteristics  to  Gas  Dynamic 
problems.  Firstly,  shock  waves  appear  frequently  in  such  problems  ana  need  to  ce 
fitted  on  an  ad  hoc  basis  into  the  characteristic  network.  Secondly,  this  network 
does  not  have  uniform  mesh  spacing,  is  non-orthogonal  ,  curvilinear  and  frequently 
highly  skewed. 

One  of  the  earliest  methods  avoiding  these  drawbacks  'was  the  Finite  Difference 
method  of  von  Neumann  and  Richtmyer  (1950).  This  applied  to  unsteady  spherical 
flow  in  Lagrangian  coordinates,  using  finite  differences  in  tne  original  ortr.ogo- 
nal  independent  variables,  and  treated  shocks  by  the  introduction  of  artificial 
viscosity.  This  was  refined  later  in  the  well-known  Lax-Wendroff ( 1954)  method. 

The  development  of  these  methods  is  well  described  in  Richtmyer  and  Morton  (1972). 


The  search  for  better  numerical  techniques  to  solve  Gas  Dynamic  problems  was  given 
further  impetus  in  the  Sputnik  era,  with  the  need  to  solve  the  flow  field  problem 
for  a  space  vehicle  on  re-entry  into  the  earth's  atmosphere,  the  so-called  Blurt 
Body  problem. 


This  led  to  a  spate  of  new  techniques,  mostly  developed  in  the  USSR.  Firstly, 
Godunov's  method  (1959)  was  presented  as  a  novel  way  to  solve  the  Lagrangian  equa¬ 
tion  in  unsteady  flow  but  grew  into  one  for  solving  the  unsteady  Culerian  equa¬ 
tions  in  two  and  three  dimensions  -yielding  the  solution  to  the  blunt  no dy  problem 
as  the  steady  flow,  asymptotic,  limit  of  the  unsteady  solution.  Godunov  used 
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discontinuity  breakdown  formulae,  in  place  of  finite  difference  formulae,  a  sig¬ 
nificant  departure  from  the  previous  approaches  of  Lax-Wendroff  and  others. 

A  solution  to  the  blunt  body  problem  predating  that  of  Godunov  was  obtained  b y 
Belotserkovski i  (1958)  who  solved  the  equations  of  steady  motion  for  symmetrical 
flow  past  a  circular  cylinder  or  sphere,  of  mixed  el  1 iptic-hyperbol ic  type,  by 
the  first  formulation  of  the  Method  of  integral  Relations  due  to  Dorodni tsyr. 
(1956).  In  this  method  moments  of  the  equations  of  motion  are  integrated  across 
the  flow  field,  the  integrands  are  represented  by  suitable  interpolation  functions 
and  the  problem  is  reduced  to  the  integration  of  ordinary  differential  equations 
in  the  coefficients  of  these  functions. 

The  steady  flow  version  of  the  Blunt  Body  problem  was  also  solved  by  Toler, in's 
method  (Gilinskii,  Telenin  and  Tinyakov,  1964),  in  which  the  unknown  is  repre¬ 
sented  by  an  interpolation  function  in  one  of  the  coordinates,  in  the  original 
equations  of  motion.  The  Method  of  Lines  (Jones  and  South,  1979)  is  a  similar 
technique  applied  over  three  or  five  mesh  points  instead  of  over  the  whole  flew 
field. 

The  last  method  developed  for  high  speed  steady  flow  problems  is  due  to  Babenko 
and  others.  The  version  applied  to  purely  supersonic  flow  is  by  Babenko, 
Voskresenski i ,  Lyubimov  and  Rusanov  (1964)  while  the  blunt  body  version  was 
developed  by  Lyubimov  and  Rusanov  (1970).  This  is  a  sophisticated  finite  dif¬ 
ference  method  applicable  to  any  steady  high  speed  inviscid  flow  problem,  especial 
ly  in  three  dimensions. 

■When  viscous  effects  are  important  we  need  to  solve  the  Navier-Stokes  equations, 
either  in  their  original,  or  approximate,  boundary  layer  form.  The  full  equations 
are  non-linear  and  elliptic.  Several  Finite  Difference  methods  have  been  develope 
for  these  but  they  become  prohibitively  expensive  as  the  governing  Reynolds  number 
is  increased.  In  the  low  speed  range  Finite  Element  methods  have  been  proposed 
for  flow  fields  limited  by  boundaries  (cavities  and  steps)  while  Spectral  methods 
have  been  used  for  investigation  of  flow  structures  in  unlimited  regions.  A  full 
account  of  these  is  given  in • Peyret-Taylor  (1983). 

In  many  applications  viscous  effects  only  need  to  be  considered  near  boundaries 
and  it  is  sufficient  to  solve  the  boundary  layer  equations  in  these  regions  in 
interaction  with  inviscid  flow  in  the  outer  regions.  Many  finite  difference  meth¬ 
ods  have  been  developed  for  the  boundary  layer  equations.  At  the  present  time, 
however,  these  can  be  replaced  by  more  recently  developed  techniques  such  as  the 
Method  of  Integral  Relations,  Finite  Element  methods,  Galerkin  techniques  and 
Spectral  Methods.  These  are  described  in  full  detail  by  Fletcher  (1983). 

The  remainder  of  this  article  deals  with  two  topics  connected  with  the  author's 
own  research.  Recent  Applications  of  the  Method  of  Integral  Relations  to  Turbulent 
and  Internal  Flows  and  Recent  Developments  in  Methods  for  Problems  in  Gas  Dynamics 
.  and  Propagation  of  Large  Amplitude  Surface  Waves. 

RECENT  APPLICATIONS  OF  THE  METHOD  OF  INTEGRAL  RELATIONS 

The  Method  of  Integral  Relations  was  first  formulated  for  Viscous  Incompressible 
Boundary  Layer  Problems  by  Dorodnitsyn  (1960).  This  formulation  was  extended  to 
compressible  flows  by  Pavlovskii  (1963)  and  to  flows  with  wall  injection  by  Liu 
(1962).  The  early  applications  were  all  to  attached  flows  but  extensions  of  the 
method  to  apply  to  separated  and  reversed  flows  were  made  by  Nielsen,  Goodwin  and 
Kuhn  (1969)  and  by  Holt  (1966,1967)  and  Holt  and  Lu  (1975). 

In  the  original  formulation  the  basic  integral  relation  is  derived  by  factoring 
the  continuity  equation  by  one  of  a  set  of  weighting  functions  f(u)  and  the 
streamwise  momentum  equation  by  its  derivative  f'(u)  where  u  is  the  streamwise 
velocity  component.  The  results  are  added  and  integrated  across  the  boundary 
layer,  using  u,  rather  than  n  (the  transverse  coordinate),  as  variable  of  integra¬ 
tion.  The  functions  f ( u)  belong  to  a  complete  set  and  taken  in  fact  as  integral 
powers  of  (1-u)  [u  is  made  dimensionless  with  respect  to  velocity  just  outside 
the  boundary  layer].  The  integrands  in  the  integral  relation  contain  the  unknown 
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transverse  velocity  gradient  3u/3n,  which,  in  the  Dorodnitsyn  formulation,  is 
represented  as  a  polynomial  in  u,  factored  by  a  term  (1-u)  to  ensure  approach 
to  zero  at  the  boundary  layer  edge.  Successive  integral  relations  in  the  nth 
approximation  yield  n  ordinary  differential  equations  for  the  coefficients  in 
the  (3u/3n)  polynomial.  The  matrix  of  these  equations  (for  the  first  deriva¬ 
tives)  is  dense  and  the  system  becomes  progressively  more  ill-conditioned  as 
the  order  of  approximation  is  increased. 

To  eliminate  this  defect  Fletcher  arid  Holt  (1975)  proposed  a  modified  formulation 
of  the  Method  of  Integral  Relations.  In  this,  the  weighting  functions  f(u) 
belong  to  a  complete  orthonormal  set  in  (0,1).  Moreover,  the  unknown  (3u/3n)  is 
expanded  in  combinations  of  the  same  orthonormal  set,  factored  by  (1-u).  The 
resulting  ordinary  differential  equations  for  the  coefficients  in  the  expansion 
then  have  a  sparse  matrix  with  diagonal  elements  only.  The  integrals  which 
occur  can  all  be  evaluated  by  quadratures  and  the  modified  version  can  be  applied 
to  attached  flows  at  any  level  of  approximation. 

The  extension  of  the  Dorodnitsyn  MIR  to  separated  flows  introduces  a  great  deal 
of  tedious  algebra,  even  in  the  lowest  order  approximation  and  it  would  be  useful 
to  extend  the  Fletcher-Hol t  modified  version  for  sue.)  flows.  The  difficulty  here 
is  that  a  factor  (u+a)'2  must  be  included  in  the  representation  of  (Gu/Gr.)  to  take 
account  of  reversed  flow  between  u  =0  (at  the  wall)  and  u  =  -a,  with  a  vertical 
tangent  in  the  ?u/3r,  -u  curve  at  u  =-a. 

A  possible  orthonormal  formulation  for  separated  flows  is  now  given,  developed 
as  a  course  term  project  at  the  University  of  California,  Berkeley  by  0.  Qzcan 
and  R.-J.  Yang.  This  uses  Chebychev  polynomials.  This  is  followed  by  an  ortho¬ 
normal  formulation  of  MIR  for  turbulent  boundary  layers,  applied  to  model  turbu- 
lent  flows  by  Yeung  and  Yang  (1981)  and  to  the  turbulent  wall  jet  problem  by 
Yang  and  Hoi t  (1983) . 


0RTH0N0RMAL  FORMULATION  OF  MIR  FOR  SEPARATING  FLOWS 


When  the  original  Method  of  Integral  Relations  is  applied  to  the  incompressible 
laminar  boundary  layer  equations  in  Dorodnitsyn  variables,  the  following  basic 
integral  relation  results 
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where  l  is  the  streamwise  coordinate,  u  the  velocity  component  in  the  ►  direction, 
0  the  reciprocal  of  (3u/3n),  n  the  transverse  coordinate  and  f|<(u)  are  weignting 
functions.  The  points  A  and  B  correspond  to  the  u  limits  of  integration  over  the 
section  of  boundary  layer  in  question.  The  values  n-|  ,02  correspond  to  n  at  A  and 
B,  respectively. 

The  velocity  field  is  divided  into  two  parts  in  the  attached  (but  retarded) 
boundary  layer  region  and  into  three  parts  in  the  separated  region. 


Attached  flow 


In  the  velocity  range  0  <u  <e,  where  e  is  slJ^lV0ic  "  0.1),  0  is  represented  by 
the  second  order  approximate  expression  ^ omrE  of  scientific  RESFJvru 
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(2) 


e  = 


(u+ot)'2(l-u) 


The  following  two  weighting  functions  are  chosen  to  simplify  the  calculation  of 
the  integrals  in  Eq.  (1) 


fk(u)  =  (c-u) ( 1-u) ( u+a) 


k+k 


,  k  = 1 ,2  . 


(3) 


Equations  (2)  and  (3)  are  substituted  in  Eq.  (1)  to  yield  two  ordinary  differen¬ 


tial  equations  for  c0  and  a 
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where  E(.  are  polynomials  in  a,  and  are  polynomials  in  a,  linear  in  c0,  with 
coefficients  containing  e. 


In  the  velocity  range  e  <u  <1,  0  is  represented  by  the  following  expression 
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(u+a)'2(  1-u)  u  -i=2  j  =  l  u ( 1  —  u )  v'u-c 
together  with  weighting  functions 

fk  *  n-u)5/4  5k  , 

where  g^  are  Chebychev  polynomials,  i  and  k  are  even,  j  is  odd. 
The  orthogonal i ty  condition  for  these  polynomials  is 
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The  limits  of  integration  e  ->-1  are  changed  to  -1  -*1  by  means  of  the  change  in 
variable 


We  can  show  that 


U  =  ~[e+l-2u]  . 
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Hence,  matching  expressions  (2)  and  (5)  at  u  =c, 


a  =  b 
o  o 


Substitution  of  Eqs.  (5)  and  (6)  in  Eq.  (1)  yields 

„  db. 

f  (N-4)  =  C’(k)  ,  k  =1 . N  , 
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where  C * ( k )  is  a  linear  combination  of  the  unknowns  bn  with  coefficients  which 
can  be  evaluated  numerically. 
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Separated  flow 


In  the  separated  flow  region  three  representations  of  3u/9n  in  terms  of  u  are 
used.  In  the  reversed  flow  range  -a  <u  <0,  close  to  the  wall,  we  use 


with  weighting  function 


a 

_ o 
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f  =  (u+a)^  . 
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In  the  intermediate  range  -a  <u  <e,  where  e  is  small  and  positive,  we  use 

a 
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Substitution  of  Eqs.  (10)  and  (11)  into  Eq.  (1)  yields  two  first  order  ordinary 
differential  equations  for  aQ  and  a. 


In  the  outer  range  c  <u  <1  we  again  use  representation  (5)  for  e  with  weighting 
functions  (6).  The  coefficients  bn  are  determined  from  the  same  equation  (9). 


This  approach  is  currently  being  applied  to  model  separated  flow  problems. 


APPLICATION  OF  MIR  TO  TURBULENT  BOUNDARY  LAYERS 


The  key  to  the  extension  of  the  Method  of  Integral  Relations  to  apply  to  turbu¬ 
lent  boundary  layer  flow  is  in  the  representation  of  the  eddy  viscosity  as  a 
function  of  mean  streamwise  velocity  component.  This  was  investigated  by  Abbott 
and  Deiwert  (1968)  and  by  Murphy  and  Rose  (1968).  Although  their' models  of  eddy 
viscosity  were  reasonable  these  were  not  well  adapted  to  the  original  formula¬ 
tion  of  MIR  and  produced  disappointing  results.  On  the  other  hand-  when  the 
modified,  orthonormal,  version  of  MIR  is  applied  to  turbulent  boundary  layers, 
with  representations  of  eddy  viscosity  similar  to  those  used  previously,  compari¬ 
sons  of  applications  to  experimental  results  in  model  cases  prove  to  be  very 
satisfactory.  This  advance  is  described  in  Yeung  and  Yang  (1981). 


The  turblent  boundary  layer  equations  in  two  dimensions  may  be  written 
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where  u  and  v  are  the  x  and  y  components  of  mean  velocity,  p  is  the  mean  pressure, 
p  the  density,  v  the  kinematic  viscosity.  The  eddy  viscosity  t  is  given  by 


-P  uV  =  c  (14) 

representing  the  Reynolds  stress. 


We  introduce  the  dimensionless  variables 
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where 
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Bernoulli's  equation  gives 
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and  Eqs.  (12)  and  (13),  in  the  dimensionless  variables,  become 
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To  apply  the  Method  of  Integral  Relations  we  introduce  a  complete  set  of  ortho¬ 
normal  functions  (over  (0,1))  f-j('J),  multiply  Eq.  (18)  by  f • ,  add  the  result  and 
integrate  with  respect  to  y  across  the  boundary  layer.  Ue  then  change  the 
variable  of  integration  from  y  to  U,  introducing  the  reciprocal  of  the  transverse 
mean  velocity  gradient 
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The  basic  integral  relation  is  then 
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The  orthonormal  functions  f^U)  are  given  by 


and 


f  (U)  =  Z  C.k(l-U)k  , 
1  k=]  1K 


1 


Vj  FU  dU 


kj 


(21) 

(22) 

(23) 


where  Sy  is  the  Kronecker  delta.  We  represent  Z  by  the  factored  orthonormal 
expansion  N_1 


Z  = 


b  +  Z  b.f.(U) 
°  .1=1  J  J 
1-U 


(24) 


The  basic  integral  relation  (21)  then  yields  the  following  set  of  ordinary  dif- 

iknc 

db. 


ferential  equations  for  the  unknown  coefficients  bQ,b^,b2,...  in  Eq.  (24): 
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These  can  be  integrated  subject  to  suitable  initial  conditions. 


The  advantages  of  the  modified  version  of  MIR  over  the  original  version  are  now 
apparent.  Firstly,  the  system  of  ordinary  differential  equations  [(25a)  and 
( 25b ) 3  can  easily  be  reduced  to  diagonal  form  while  the  corresponding  equations 
in  the  original  formulation  are  highly  coupled.  Secondly,  the  representation 
(24)  for  Z  allows  for  greater  flexibility  in  reproducing  the  highly  inflected 
velocity  profile  characteri zing  turbulent  (as  opposed  to  laminar)  boundary 
layers.  Thirdly,  in  the  orthonormal  version,  the  integral 

,1 

(1  + 
o 


-) 

u 


fV(U)dU 


can  be  evaluated  by  quadrature,  while  in  the  original  version  this  must  be 
reduced  to  an  algebraic  expression.  This  is  a  tedious  task,  increasingly  dif¬ 
ficult  to  carry  out  as  the  order  of  approximation  is  increased,  since  the  term 


c/p  is  represented  by  a  complicated  combination  of  exponentials  and  powers  in  U. 


Turbulence  modelling 


The  eddy  viscosity  c/p  is  represented  by  two  formulae,  one  based  on  a  model  due 
to  Spalding  (1961)  and  Kleinstein  (1967)  applicable  near  the  wall  and  the  other, 
applicable  in  the  outer  part  of  the  boundary  layer,  based  on  the  wake  model  of 
Clauser  (1956).  The  value  of  U  at  the  junction  of  these  two  representations  is 
denoted  by  Um  and  is  determined  by  conditions  that  the  e/r  -U  curve  should  be 
continuous  with  continuous  slope  at  U  =Um. 


In  terms  of  U  and  Z  the  eddy  viscosity  is  represented  by: 
For  0  <  U  <U 
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~  =  0.04432[e 


0.4U  >/U„Z„Re1/2 


e  o 


-1  -0.4U  ’'u  Z  Re1/2  -0.08U2U  Re1/2Zj  . 
e  o  e  o 


For  U  <  U  <  1 
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-  =  0.0168  U  Re1/2  (  (l-U)ZdU  . 

w  e  J0 
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(27) 


The  value  of  Um  is  determined  from 


0.04432[e 


0.4U  MjTzRe1/2 


e  o 


-1 


0.4U  4z  Re1/2  =  0.08U2U  Z  Re1/2 
e  o 


e  o 


=  0.0168  UgRe1/2  j1  (l-U)ZdU 


(28) 


Yeung  and  Yang  (1981)  applied  this  formulation  to  three  model  flows  based  on  data 
presented  at  the  1968  Stanford  Conference  on  Turbulent  Flows  (Coles  and  Hirst 
(1968)).  The  first  corresponds  to  a  zero  pressure  qradient  (identified  as  ID 
1400  in  the  Stanford  proceedings),  the  second  to  adverse  pressure  gradient  (ID 
1100)  and  the  third  to  favorable  pressure  gradient  (ID  1300).  In  the  first  case 
good  results  were  obtained  with  the  third  approximation  N  =3  while  agreement  with 
experiment  was  excellent  for  N  =4,  N  =5.  The  CPU  time  for  N  =4  was  only  8  secs 
(usinq  a  CDC  7600  computer).  In  the  third  case  sufficient  accuracy  was  obtained 
with  N=3  and  higher  order  approximations  were  not  needed.  For  adverse  pressure 
gradient  flow  results  were  partially  satifactory  for  N  =5  but  it  is  evident  that 


7. 
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the  e/p  model  needs  to  be  improved  in  this  area. 
The  wal 1  jet  probl em 


Yang  and  Holt  (1983)  used  the  orthonormal  version  of  MIR  for  turbulent  boundary 
layers  to  investigate  the  effect  of  injecting  a  parallel  stream  on  turbulent 
flow  in  a  pipe.  This  has  important  applications  to  coal  gasification  systems, 
in  which  certain  products  (especially  sulfides  and  oxides)  can  cause  serious  cor¬ 
rosion  at  the  walls  of  pipes  in  the  heat  exchange  section. 

To  simplify  the  problem  the  flow  field  is  divided  into  three  parts.  Firstly, 
just  downstream  of  the  station  where  the  wall  jet  is  introduced,  turbulent  mixing 
between  the  uniform  flow  in  the  pipe  and  uniform  flow  in  the  parallel  jet  is 
investigated.  The  injected  gas  and  main  stream  gas  contain  different  species  and 
both  turbulent  and  molecular  diffusion  are  accounted  for.  Below  the  mixing 
region  a  turbulent  boundary  layer  develops  along  the  wall  of  parallel  jet  and  its 
growth  is  determined  by  a  separate  application  of  MIR.  This  second  cart  of  the 
flow  field  is  not  initially  affected  by  the  nixing  process.  The  third  part  of 
the  flow  field  extends  downstream  of  the  station  where  the  wall  boundary  layer 
and  lower  part  of  the  mixing  layer  intersect.  The  velocity  profile  in  this 
region  of  interaction  initially  has  two  inflexion  points  but  the  orthonorr.al 
version  of  MIR  is  sufficiently  flexible  to  permit  a  faithful  representation  of 
the  Z-U  behavior,  and  the  growth  of  this  interacting  flow  until  a  full  turbulent 
boundary  layer  is  developed  and  calculated. 

The  calculation  establishes  that  introduction  of  a  wall  jet  can  protect  a  pine 
wall  from  corrosive  effects  for  up  to  100  jet  thicknesses  downstream. 

NUMERICAL  METHODS  FOR  UNSTEADY  FLOW  PROBLEM, S 

The  equations  of  motion  of  one  dimensional  unsteady  flow  are  easily  reduced  to 
characteristic  form  and  the  earliest  numerical  methods  for  solving  problems 
governed  by  such  equations  were  based  on  the  Method  of  Characteristics.  The 
latter  has  the  disadvantage  of  being  tied  to  a  non-orthogonal  coordinate  network 
which  has  to  be  built  up,  step  by  step,  as  the  numerical  procedure  is  advanced. 

To  overcome  this  drawback  Godunov  (1959)  proposed  a  numerical  scheme  for  solving 
Lagrangian  equations  of  motion  which  could  be  executed  in  the  physical  x  (dis¬ 
tance),  t  (tine)  plane.  He  later  extended  this  scheme  to  apply  to  the  Eulerian 
equations  in  one  or  more  dimensions  and  used  this  to  solve  the  Blunt  Body 
Problem.  These  early  Godunov  schemes  are  of  first  order  accuracy  and  are 
monotonic  in  character-,  that  is,  if  the  initial  form  of  the  unknown  is  mono- 
tonically  increasing  this  property  is  preserved  as  the  Godunov  schemes  are 
applied  at  successive  time  intervals.  Godunov  (1970)  subsequently  introduced  a 
scheme  of  second  order  accuracy,  based  on  a  predictor-corrector  approach.  This 
is  not  monotonic  in  character  but  has  been  successfully  applied  to  many  problems 
in  Gas  Dynamics  with  plane  symmetry  and  in  Shallow  Water  Theory.  In  the  period 
between  the  publication  of  the  two  Godunov  schemes,  Glimm  (1965)  proposed  a 
modification  of  the  first  Godunov  scheme  which,  in  principle,  has  second  order 
accuracy.  We  now  give  a  brief  comparison  between  Godunov's  first  scheme  and 
Glimm's  scheme,  using  explanations  of  the  latter  by  Chorin  (1976,1977).  The 
comparison  is  made  with  reference  to  the  one  dimensional  acoustic  equations  and 
is  a  summary  of  the  discussion  in  Holt  (1984). 

Solution  of  the  acoustic  wave  equation 


Consider  the  one  dimensional  acoustic  wave  equation 

3U  +  ,  3M  -  o 

9t  9x  ° 
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where  a  is  a  positive  constant  (the  acoustic  speed)  with  initial  conditions 


U  =  f(x)  ,  t  =  0  .  (30) 

The  analytical  solution  to  this  equation  is 

U  =  f(x-at)  .  (31) 

If  the  initial  wave  form  is  a  Heaviside  unit  function 


the  general  solution  is 


f(x) 

=  0 

,  x  <  0 

f  (  x) 

=  1 

,  x  >  0 

U(x,t) 

=  0 

,  x  <  at 

U(x,t) 

=  1 

,  x  >  at 

(32) 


(33) 


representing  a  step  function  propagated  unchanged  along  the  characteristic  line 
x  =at. 


To  solve  problem  (29),  (30)  by  Godunov's  first  scheme  we  divide 
a  series  of  small  cells  (usually  of  equal  spacinq)  and  represen 
staircase  function  so  that  f(x)  is  constant  in  each  cell.  At  e 
there  is  a  discontinuity  in  the  initial  value  of  U  and,  for  t  > 
problem  of  breakdown  in  discontinuity.  In  other  words,  we  solv. 
problems  (23),  (32).  The  same  procedure  is  followed  at  all  late. 


e  x  axis  into 
'x)  by  a 
cell  bouncary 
we  solve  the 
a  accession  of 
Line  intervals. 


Thus,  at  a  general  time  (after  n  time  steps  each  of  duration  k)  we  solve  the 
initial  value  problem:  Solve  (29)  with 


U  =  0  ih  <x  <(i  +'j)h  t  =nk 

U  =  1  (i  +!i)h  <  x  <  ( i  + 1  )h  t  =nk 


(34) 


If  we  transfer  the  origin  {(i  +'j)h,  nk)  to  (0,0)  this  problem  has  the  solution 

U(x,t)  =0  X  <  aT 

( 

U(x,t)  -  1  X  >  aT  , 


where  X,T  are  coordinates.  In  the  Godunov  scheme  we  always  use  this  solution  at 
the  cell  boundary  itself.  To  satisfy  stability  of  the  Godunov  scheme  we  require 

ah  >  k  (36) 

so  that  the  boundary  {(i  +l)h,  (n  +!*)k}  is  to  the  right  of  X  -aT  =0.  Therefore, 

U  +1  =  1  .  (37) 

n  +*s 

so  that  the  Godunov  scheme  always  moves  the  step  jump  0-1  a  distance  h  to  the 
right.  As  a  consequence  it  can  be  shown  that,  after  N  whole  time  steps  the 
Godunov  scheme  causes  the  step  jump  to  have  moved  a  distance 

{h/k  -a}T  (38) 

beyond  its  position  as  given  by  the  analytical  solution. 


In  Glimm's  scheme  we  record  the  solution  to  tne  breakdown  problem  (29), (34)  over 
a  half  step  range  or.  either  side  of  the  cell  boundary.  We  then  sample  this  solu¬ 
tion  at  a  randomly  chosen  point  within  the  range  (-h/2, h/2)  in  X.  The  stability 
(Courant-Friedrichs-Lewy)  condition  ensures  that  the  characteristic  line  X  =a^ 
intersects  the  line  T  =hk  inside  (-h/2,h/2).  The  greater  flexibility  in  the 
G1  imm  scheme,  as  compared  with  the  first  Godunov  scheme,  permits  us  to  sample  a 
certain  number  of  intersection  points  in  (-hh,ph)  on  both  sides  of  X=aT.  If 
the  randomly  chosen  samples  are  uniformly  distributed  over  ( -h/2 , h/2 )  it  car  bc- 
shown  that  balance  between  samples  to  the  left  and  those  to  the  right  of  X  =aT  is 
that  required  to  ensure  that  the  path  of  tne  discontinuity  0-1  stays  on  course. 
In  this  sense,  then,  the  Glimm  scheme  is  more  accurate  than  the  first  Godunov 
scheme. 

Shallow  water  wave  propagation 

Li  and  Holt  (1981)  applied  Glimm's  method  to  the  problem  of  shallow  water  waves 
generated  by  large,  near  surface,  disturbances.  The  calculation  is  an  extension 
of  the  work  by  Sod  (1977)  on  spherical,  or  cylindrical,  explosions  in  gases  and 
the  paper  by  Flores  and  Holt  (1981)  on  underwater  explosions.  Li  and  Holt  first 
tested  Glimm's  method  on  the  classical  Dam  Break  Problem  and  then  applied  it  to 
the  Dam  Break  Problem  with  cylindrical  symmetry.  It  was  then  used  to  calculate 
large  amplitude  surface  waves  generated  by  near  surface  explosions. 

The  shallow  water  equations,  omitting  dissipative  terms,  may  be  written 

Ut  +{F(U)lr  =  -W(U)  (39) 

where 

i(u/r)(n  +d) 

W(U)  = 

0 

and 

i  =  0  for  plane  symmetry 

i  =  1  for  cylindrical  symmetry 

In  Eqs.  (40)  d  is  the  undisturbed  ocean  depth,  n  is  the  displacement  of  the  ocean 
surface  from  its  undisturbed  position,  r  is  the  space  coordinate. 

Glimm's  method  is  only  applicable  to  equations  of  motion  in  conservation  form. 
Following  Sod  (1977)  we  therefore  solve  Eqs.  (39)  by  a  splitting  method.  We  soiv 
the  homogeneous  equations 

Ut+{F(U))r  =  0  (41) 

by  the  established  Glimm  technique  for  plane  flow  equations  and  subsequently 
determine  the  non-homogeneous  term  from 

Ut  »  -W(U)  .  (42) 

To  apply  Glimm's  method  to  Eqs.  (41)  we  use  h  as  a  constant  space  interval  and  k 
as  the  constant  time  interval.  After  time  t  =nk  we  represent  the  solution  by  a 
staircase  function  u(?  such  that 

U(r,nk)  =  u"+]  r  >  (i  +Jj)h  , 

U(r,nk)  =  u"  r  <  (i  +  ij)b  . 

J 
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(A3) 


Thus  U  is  constant  in  each  cell,  with  spacing  h,  along  the  r  axis  and  jur.os  dis- 
continuously  across  each  cell  boundary.  To  determine  U  at  t  =(n  Jk  we  allow 
each  cell  boundary  to  be  suddenly  removed  at  t  rnk  and  solve  a  series  of  Riemann 
breakdown  problems  to  determine  the  solution  in  nk  <  t  <  (n  +'Jk.  In  Glimm's 
method  we  sample  this  solution  at  a  randomly  chosen  point  in  the  ranae  (-;>.,'.h) 
on  either  side  of  a  cell  boundary.  If  0  is  an  equidistributed  random  variable  in 
(~Wj)  Glimm's  method  then  gives 

u"+!j2  =  U{(i  +6)h,  (n  +>i)k}  (44) 

The  grid  and  the  wave  fronts  used  in  solving  successive  Riemann  problems  are 
shown  in  Fig.  1,  while  the  sampling  procedure  is  shown  in  Fig.  2. 

After  solving  the  Riemann  problem  at  t  -(n  +  '2)k  the  solution  (43)  is  suostituted 
in  the  right-hand  side  of  Eq.  (42)  which  is  solved  as  a  simple  difference  equation 

n-f1, 

for  a  corrected  u-,,*  .  The  corrected  values  are  used  as  initial  data  to  apply  the 
1+2 

Glimin  method  in  the  next  half  time  interval  (n  +  R)k  <  t  <  nk.  The  Glimm  solution 
at  t  =nk  is  again  corrected  from  Eq.  (42). 

The  solution  of  the  Riemann  problem  at  each  cell  boundary  is  given  by  algebraic 
formulae,  the  details  of  which  depend  on  the  nature  of  the  discontinuities  across 
the  boundary.  Waves  are  propaciated  in  botii  directions,  when  each  discontinuity 
breaks  down,  and  can  be  either  expansion  waves  or  bores  (the  shallow  water  equiva¬ 
lent  of  shocks).  The  Riemann  problem  solutions  are  given  in  full  in  Li  and  Holt 
( 198h) . 

The  initial  conditions  for  the  classical  Dam  Break  Problem  are  shewn  in  Fig.  3 
(Fig.  4, Li  and  Holt).  The  dam  maintains  the  drop  in  water  level  shown.  It  is 
suddenly  removed  at  time  t  =0.  Subsequently  a  bore  is  propagated  to  the  left 
while  an  expansion  is  propagated  to  the  left.  In  tne  Holt-Li  calculations  a 
solid  wall  is  introduced  to  ttie  left  of  the  dam  so  that  reflection  of  the  bore  can 
be  calculated.  The  solutions  before  and  after  reflection  are  shown  in  Figs.  4  arc 
4.  These  calculations  were  made  by  the  Glimin  method  alone,  since  Eq.  (42)  was  net 
required  in  the  plane  flow  case.  The  Holt-Li  method  was  then  applied  to  the 
cylindrical  dam  break  problem.  In  this  case  a  cylindrical  bore  converged  cn  its 
center  and  was  reflected  there,  in  a  manner  similar  to  a  cylindrical  implosion  in 
a  gas.  The  wave  profiles  in  this  case  are  shown  in  Fig.  6. 

The  final  calculation  deals  with  the  Upper  Critical  Depth  problem.  This  concerns 
the  generation  of  large  amplitude  surface  waves  by  spherical  explosions  detonated 
at  different  depths  below,  but  near  to,  tne  ocean  surface.  The  pressure  field 
from  the  underwater  explosion  was  calculated  previously  by  Falade  and  Holt  (1373). 
The  data  were  used  to  provide  initial  conditions  in  the  surface  wave  calculation 
using  Glimm's  method.  Figure  7  shows  the  most  important  result  of  this  calcula¬ 
tion  namely,  that  for  different  depths  of  explosive  charge  the  maximum  amplitude 
of  surface  wave  generated  occurs  when  the  charge  depth  is  equal  to  one  half  of 
the  charge  radius.  This  confirms  the  Upper  Critical  Depth  phenomenon  observed  in 
field  tests. 
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Fig.  1.  Sequence  of  Riemann  problems  on  qrid. 
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- EXACT  SOLUTION 

A  NUMERICAL  SOLUTION 
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Fig.  4.  Dam  break  problem  with  plane  symmetry. 
Time  t  =0.69. 
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-  EXACT  SOLUTION 

A  NUMERICAL  SOLUTION 
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g.  5.  Dam  break  problem  with  plane  symmetry. 
T  i  me  t  =  1 .  36 . 
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Dam  break  problem  with  cylindrical  symmetry. 
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